Здесь самостоятельно реализован алгоритм кластеризации KMeans.
Алгоритм K-Means.
Параметры:
Атрибуты:
Метод predict(): Новый объект определяется в кластер, центр которого расположен ближе всех к этому объекту.
from sklearn.base import BaseEstimator, ClusterMixin
import math
#from sklearn.cluster import _kmeans_plusplus
class My_Clust(BaseEstimator, ClusterMixin):
"""
My implementation of clustering algorithm
"""
def __init__(self, k=8, init='random', n_init=10, max_iter=300, tol=0.0001):
"""
Please add necessary algoritm parameters to class constructor.
"""
self.k = k
self.init = init
self.n_init = n_init
self.max_iter = max_iter
self.tol = tol
self.inertia_ = 0.0
self.cluster_centers_ = [[]]
self.labels_ = []
def _find_init_cluster_centers(self, X, init):
n_samples = len(X)
n_clusters = self.k
if isinstance(init, str) and init == 'random':
centers = X[np.random.permutation(n_samples)][:n_clusters]
elif isinstance(init, str) and init == 'kmeans++':
centers, _ = kmeans_plusplus(X, n_clusters)
return centers
def _make_labels(self, X, centers):
mapping = []
features = len(X[0])
inertia = 0
for s in X:
distances = []
inertia_list = []
for c in centers:
distance = sum([(s[p] - c[1][p]) ** 2 for p in range(features)])
distances.append(np.sqrt(distance))
inertia_list.append(distance)
min_distance = distances.index(min(distances))
mapping.append((s, min_distance))
min_square_distance = inertia_list.index(min(inertia_list))
inertia += min_distance
return mapping, inertia
def _compute_centers(self, labels):
n_clusters = self.k
new_centers = []
n_features = len(labels[0][0])
for n in range(n_clusters):
cluster = list(filter(lambda x: x[1] == n, labels))
new_centers.append((n, np.array([np.mean([x[0][p] for x in cluster]) for p in range(n_features)])))
return new_centers
def fit(self, X):
"""
Use data matrix X to compute model parameters
"""
self.n_init += 1
self.cluster_centers_ = self._find_init_cluster_centers(X, self.init)
centers = [(i, c) for i, c in enumerate(self.cluster_centers_.copy())]
for _ in range(self.max_iter):
labels, inertia = self._make_labels(X, centers)
new_centers = self._compute_centers(labels)
center_distanes = np.sqrt(sum([(centers[p][1] - new_centers[p][1]) ** 2 for p in range(len(centers))]))
if all([d <= self.tol for d in center_distanes]):
self.cluster_centers_ = np.array([c[1] for c in centers])
self.labels_ = [l[1] for l in labels]
self.inertia_ = inertia
return
centers = new_centers.copy()
def predict(self, X):
"""
Using computed model parameters predict cluster
for all objects from matrix X
return
"""
centers = [(i, c) for i, c in enumerate(self.cluster_centers_)]
labels_, self.inertia_ = self._make_labels(X, centers)
self.labels_ = np.array([x[1] for x in labels])
return self.labels_
Результаты кластеризации должны совпадать кастомного и библиотечного алгоритмов должны совпадать.
import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.cluster import KMeans
iris = load_iris()
X = iris.data # использовать для кластеризации
y = iris.target # истинные метки цветков
kmeans = KMeans(n_clusters=3, init='random')
clust = My_Clust(k=3)
kmeans.fit(X)
clust.fit(X)
kmeans.cluster_centers_, clust.cluster_centers_
(array([[6.85 , 3.07368421, 5.74210526, 2.07105263],
[5.9016129 , 2.7483871 , 4.39354839, 1.43387097],
[5.006 , 3.428 , 1.462 , 0.246 ]]),
array([[6.85384615, 3.07692308, 5.71538462, 2.05384615],
[5.006 , 3.428 , 1.462 , 0.246 ],
[5.88360656, 2.74098361, 4.38852459, 1.43442623]]))
X = pd.DataFrame(X, columns=iris.feature_names)
X['class'] = [iris.target_names[i] for i in y]
sns.pairplot(X, hue='class', plot_kws={'alpha':0.5}, vars=iris.feature_names)
plt.show()
X['class'] = kmeans.labels_
sns.pairplot(X, hue='class', plot_kws={'alpha':0.5}, palette="tab10")
plt.show()
X['class'] = clust.labels_
sns.pairplot(X, hue='class', plot_kws={'alpha':0.5}, palette="tab10")
plt.show()
Анализ датасета по различным городам США.
pd.set_option('display.max_colwidth', None)
data_desc = pd.read_csv('https://raw.githubusercontent.com/katearb/files/main/Data_usa_city_Description.txt', sep=':')
data_desc
| Attribute | Description | |
|---|---|---|
| 0 | Place | City, state (postal code) |
| 1 | Climate & Terrain | Very hot and very cold months, seasonal temperature variation, heating- and cooling-degree days, freezing days, zero-degree days, ninety-degree days. |
| 2 | Housing | Utility bills, property taxes, mortgage payments. |
| 3 | Health Care & Environment | Per capita physicians, teaching hospitals, medical schools, cardiac rehabilitation centers, comprehensive cancer treatment centers, hospices, insurance/hospitalization costs index, flouridation of drinking water, air pollution. |
| 4 | Crime | Violent crime rate, property crime rate. |
| 5 | Transportation | Daily commute, public transportation, Interstate highways, air service, passenger rail service. |
| 6 | Education | Pupil/teacher ratio in the public K-12 system, effort index in K-12, accademic options in higher education. |
| 7 | The Arts | Museums, fine arts and public radio stations, public television stations, universities offering a degree or degrees in the arts, symphony orchestras, theatres, opera companies, dance companies, public libraries. |
| 8 | Recreation | Good restaurants, public golf courses, certified lanes for tenpin bowling, movie theatres, zoos, aquariums, family theme parks, sanctioned automobile race tracks, pari-mutuel betting attractions, major- and minor- league professional sports teams, NCAA Division I football and basketball teams, miles of ocean or Great Lakes coastline, inland water, national forests, national parks, or national wildlife refuges, Consolidated Metropolitan Statistical Area access. |
| 9 | Economics | Average household income adjusted for taxes and living costs, income growth, job growth. |
| 10 | Longitude | Longitude |
| 11 | Latitude | Latitude |
| 12 | Population | Population |
Housing и Crime - наоборот.Population- статистический признак, не имеющий интерпретации как “лучше-хуже”.Place - уникальный идентификатор объекта (города), он не должен использоваться при кластеризации.Longitude и Latitude. Их также не следует использовать при кластеризации данных.data = pd.read_csv('https://raw.githubusercontent.com/katearb/files/main/Data_usa_city.txt', sep=' ')
data
| Place | Climate | HousingCost | HlthCare | Crime | Transp | Educ | Arts | Recreat | Econ | Long | Lat | Pop | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Abilene,TX | 521 | 6200 | 237 | 923 | 4031 | 2757 | 996 | 1405 | 7633 | -99.6890 | 32.5590 | 110932 |
| 1 | Akron,OH | 575 | 8138 | 1656 | 886 | 4883 | 2438 | 5564 | 2632 | 4350 | -81.5180 | 41.0850 | 660328 |
| 2 | Albany,GA | 468 | 7339 | 618 | 970 | 2531 | 2560 | 237 | 859 | 5250 | -84.1580 | 31.5750 | 112402 |
| 3 | Albany-Schenectady-Troy,NY | 476 | 7908 | 1431 | 610 | 6883 | 3399 | 4655 | 1617 | 5864 | -73.7983 | 42.7327 | 835880 |
| 4 | Albuquerque,NM | 659 | 8393 | 1853 | 1483 | 6558 | 3026 | 4496 | 2612 | 5727 | -106.6500 | 35.0830 | 419700 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 324 | Worcester,MA | 562 | 8715 | 1805 | 680 | 3643 | 3299 | 1784 | 910 | 5040 | -71.7950 | 42.2720 | 402918 |
| 325 | Yakima,WA | 535 | 6440 | 317 | 1106 | 3731 | 2491 | 996 | 2140 | 4986 | -120.5130 | 46.5950 | 172508 |
| 326 | York,PA | 540 | 8371 | 713 | 440 | 2267 | 2903 | 1022 | 842 | 4946 | -76.7280 | 39.9600 | 381255 |
| 327 | Youngstown-Warren,OH | 570 | 7021 | 1097 | 938 | 3374 | 2920 | 2797 | 1327 | 3894 | -80.7290 | 41.1700 | 531350 |
| 328 | Yuba-City,CA | 608 | 7875 | 212 | 1179 | 2768 | 2387 | 122 | 918 | 4694 | -121.6220 | 39.1280 | 101979 |
329 rows × 13 columns
Выполнить необходимую предобработку данных.
Выполнить кластеризацию иерархическим методом. (агломеративно)
Выполнить кластеризацию методом Dbscan.
Выполнить кластеризацию методом kmeans.
Выполнить кластеризацию другими методами.
Сравнить результаты.
Выбрать одно разбиение, наиболее подходящее на ваш взгляд. Проинтерпритировать.
Оценить, как полученные кластеры распределены географически.
!pip install folium
Collecting folium Using cached folium-0.12.1-py2.py3-none-any.whl (94 kB) Requirement already satisfied: numpy in c:\users\public\anaconda\lib\site-packages (from folium) (1.20.1) Requirement already satisfied: jinja2>=2.9 in c:\users\public\anaconda\lib\site-packages (from folium) (2.11.3) Requirement already satisfied: requests in c:\users\public\anaconda\lib\site-packages (from folium) (2.25.1) Collecting branca>=0.3.0 Using cached branca-0.4.2-py3-none-any.whl (24 kB) Requirement already satisfied: MarkupSafe>=0.23 in c:\users\public\anaconda\lib\site-packages (from jinja2>=2.9->folium) (1.1.1) Requirement already satisfied: chardet<5,>=3.0.2 in c:\users\public\anaconda\lib\site-packages (from requests->folium) (4.0.0) Requirement already satisfied: idna<3,>=2.5 in c:\users\public\anaconda\lib\site-packages (from requests->folium) (2.10) Requirement already satisfied: certifi>=2017.4.17 in c:\users\public\anaconda\lib\site-packages (from requests->folium) (2020.12.5) Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\users\public\anaconda\lib\site-packages (from requests->folium) (1.26.4) Installing collected packages: branca, folium Successfully installed branca-0.4.2 folium-0.12.1
from sklearn.preprocessing import MinMaxScaler
from sklearn.cluster import AgglomerativeClustering, DBSCAN
from sklearn.metrics import silhouette_score
from sklearn.mixture import GaussianMixture
from sklearn.neighbors import NearestNeighbors
from scipy.cluster.hierarchy import fcluster, linkage, dendrogram
import folium
df = data.drop(['Place', 'Long', 'Lat'], axis=1)
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 329 entries, 0 to 328 Data columns (total 10 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Climate 329 non-null int64 1 HousingCost 329 non-null int64 2 HlthCare 329 non-null int64 3 Crime 329 non-null int64 4 Transp 329 non-null int64 5 Educ 329 non-null int64 6 Arts 329 non-null int64 7 Recreat 329 non-null int64 8 Econ 329 non-null int64 9 Pop 329 non-null int64 dtypes: int64(10) memory usage: 25.8 KB
scaler = MinMaxScaler()
df = pd.DataFrame(scaler.fit_transform(df), columns=df.columns)
for m in ['euclidean', 'minkowski', 'cityblock', 'seuclidean']:
for l in ['weighted', 'single']:
print(m, l)
distance_matrix = linkage(df, metric=m, method=l)
dn = dendrogram(distance_matrix)
plt.show()
for l in ['ward', 'centroid']:
print('euclidean', l)
distance_matrix = linkage(df, metric='euclidean', method=l)
dn = dendrogram(distance_matrix)
plt.show()
euclidean weighted
euclidean single
minkowski weighted
minkowski single
cityblock weighted
cityblock single
seuclidean weighted
seuclidean single
euclidean ward
euclidean centroid
кажется, что Евклидово расстояние в совокупности с методом ward дает самый разумный результат из всех. Визуально четко выделяются 3-4 кластера. В то время как в остальных случаях, дендрограммы либо не способны выделить какие-то ни было кластеры, либо очень сильно их искажают.
K = range(2, 8)
scores = []
for k in K:
agglomerative = AgglomerativeClustering(linkage='average', n_clusters=k).fit(df)
labels = agglomerative.labels_
scores.append(silhouette_score(df, labels, metric = 'euclidean'))
print(scores)
plt.plot(K, scores)
plt.title('silhouette score')
plt.show()
[0.6342602660378628, 0.4096425727541542, 0.4066173623108926, 0.3327832822547907, 0.31280265675863367, 0.28100792432687166]
Оптимальное количество кластеров = 3
agglomerative = AgglomerativeClustering(n_clusters=3, affinity='euclidean')
agglomerative.fit(df)
AgglomerativeClustering(n_clusters=3)
neigh = NearestNeighbors(n_neighbors=5)
nbrs = neigh.fit(df)
distances, _ = nbrs.kneighbors(df)
distances = np.sort(distances, axis=0)
distances = distances[:,1]
plt.plot(distances)
plt.xlabel('sorted distances')
plt.ylabel('mean distance over 3 near neighbours')
plt.show()
Оптимальное eps = 0.4
dbscan = DBSCAN(eps=0.4, min_samples=5, leaf_size=60)
dbscan.fit(df)
DBSCAN(eps=0.4, leaf_size=60)
kmean = KMeans(n_clusters=3)
kmean.fit(df)
KMeans(n_clusters=3)
gmm = GaussianMixture(n_components=3, init_params='kmeans')
gmm.fit(df)
GaussianMixture(n_components=3)
agg_sil = silhouette_score(df, agglomerative.labels_, metric = 'euclidean')
print('silhouette score', agg_sil)
silhouette score 0.13781545274723117
print('silhouette score = 1, так как у нас всего 1 кластер')
silhouette score = 1, так как у нас всего 1 кластер
km_sil = kmean.labels_
silhouette_score(df, km_sil, metric = 'euclidean')
0.17996676101155756
gmm_sil = gmm.predict(df)
silhouette_score(df, gmm_sil, metric='euclidean')
0.0667132133661567
По значениям Силуэта видно, что все кластеризации примерно одинаково спраились с задачей (за исключением dbscan, который выделил только 1 кластер и выбросы). Лучше всех справился Gaussian Mixture.
|
|
Оценим распредление точек по кластерам с помощью pairplot.
df['class'] = agglomerative.labels_
sns.pairplot(df, hue='class', plot_kws={'alpha':0.5}, palette="tab10")
plt.show()
df.drop(['class'], axis=1, inplace=True)
На большинстве признаков кластеры перемешаны между собой. Наиболее отчетливо кластеры различаются в строке TRansp и Pop: 0 кластер всегда имеет более маленькие значения, по сравнению с кластером 2. Кластер 1 больше похож на выбросы, поскольку точки внутри этого кластера сильно разбросаны, плотность мала. Также хорошо видно расределение на паре признаков Crime-HousingCost.
df['class'] = dbscan.labels_
sns.pairplot(df, hue='class', plot_kws={'alpha':0.5}, palette="tab10")
plt.show()
df.drop(['class'], axis=1, inplace=True)
DBSCAN смог только выделить шум.
df['class'] = kmean.labels_
sns.pairplot(df, hue='class', plot_kws={'alpha':0.5}, palette="tab10")
plt.show()
df.drop(['class'], axis=1, inplace=True)
Кластеризация kmeans очень похожа на агломеративную, за исключением, что 1 кластер вбирает в себя более количество точек, что визуально выглядит более правильно и логично.
df['class'] = gmm.predict(df)
sns.pairplot(df, hue='class', plot_kws={'alpha':0.5}, palette='tab10')
plt.show()
Результаты Gаussian Mixture схожи с результатами kmeans, однако немного больше точек, определенные kmeans в 1 кластер, отнесены здесь ко 2. Из-за чего значение силуэта немного ниже.
Результаты Gaussian Mixture кажутся мне наиболее правильными на этом наборе данных.
Точки 0 кластера часто располагаются в левом нижнем углу графика, то есть значения точек для определенных признаков в среднем меньше, чем значения точек двух других кластеров. В это же время, 0 и 2 кластеры на некоторых признаках никак дистибутируются: их значения накладываются друг на друга. Это привело к тому, что DBSCAN объединил всё в 1 кластер.
2 кластер чаще всего находит между 1 и 0 кластерами, то есть его элементы обладают средними значениями по признакам. Это самый много численный класс.
1 кластер -- очень разряжен, его элементы обладают, как правило, наибольшими значениями по сранению с остальными кластерами (За исключением CRime и Econ). Этот кластер самый малочисленный. Некоторые его элементы были определены дбсканом как выбросы, что логично, ведь на графике они сильно удалены от общей массы точек. Агломеративаная кластеризация его многие точки отнесла к 2 кластеру.
Таким образом:
2 кластер -- города со средней численностью населения, средним (не жарким, не холодным) климатом, уровнем здравоохранения, часто хорошими возмодностями образования, уровнем культурного рейтинга, однако довольно высоким уровнем преступлений и малым количеством рекреационных объектов.
0 кластерм -- города с маленькой численностью населения, низкими затратами на жкх и тп., низким уровнем развития здравоохранения, низкими уровнем преступности, слаборазвитыми системами образования, культуры, рекреационных возможностей.
1 кластер -- крупные города, с большими возможностями для культурного развития и развлечений, высокими ценами жкх, развитой системой здравоохранения, средним уровнем преступности, часто с хорошими возможностями для образования.
df['Place'] = data['Place'].apply(lambda x: x.split(',')[-1].split('-')[0])
df
| Climate | HousingCost | HlthCare | Crime | Transp | Educ | Arts | Recreat | Econ | Pop | class | Place | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.516770 | 0.056328 | 0.024849 | 0.280822 | 0.385829 | 0.507692 | 0.016651 | 0.245556 | 0.661572 | 0.005859 | 0 | TX |
| 1 | 0.583851 | 0.161193 | 0.206609 | 0.263927 | 0.499733 | 0.354327 | 0.097225 | 0.518222 | 0.188176 | 0.072759 | 1 | OH |
| 2 | 0.450932 | 0.117959 | 0.073652 | 0.302283 | 0.185294 | 0.412981 | 0.003263 | 0.124222 | 0.317952 | 0.006038 | 0 | GA |
| 3 | 0.460870 | 0.148747 | 0.177789 | 0.137900 | 0.767112 | 0.816346 | 0.081192 | 0.292667 | 0.406489 | 0.094136 | 1 | NY |
| 4 | 0.688199 | 0.174991 | 0.231843 | 0.536530 | 0.723663 | 0.637019 | 0.078387 | 0.513778 | 0.386734 | 0.043458 | 1 | NM |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 324 | 0.567702 | 0.192414 | 0.225695 | 0.169863 | 0.333957 | 0.768269 | 0.030551 | 0.135556 | 0.287671 | 0.041414 | 1 | MA |
| 325 | 0.534161 | 0.069314 | 0.035097 | 0.364384 | 0.345722 | 0.379808 | 0.016651 | 0.408889 | 0.279885 | 0.013357 | 0 | WA |
| 326 | 0.540373 | 0.173800 | 0.085820 | 0.060274 | 0.150000 | 0.577885 | 0.017110 | 0.120444 | 0.274117 | 0.038776 | 1 | PA |
| 327 | 0.577640 | 0.100752 | 0.135007 | 0.287671 | 0.297995 | 0.586058 | 0.048419 | 0.228222 | 0.122422 | 0.057053 | 1 | OH |
| 328 | 0.624845 | 0.146962 | 0.021647 | 0.397717 | 0.216979 | 0.329808 | 0.001235 | 0.137333 | 0.237779 | 0.004768 | 0 | CA |
329 rows × 12 columns
pd.crosstab(df['Place'], df['class']).plot(kind='pie', subplots=True, figsize=(20, 20), legend=False)
plt.show()
pd.crosstab(df['Place'], df['class']).plot(kind='bar', figsize=(30, 8), legend=True, stacked=True)
<AxesSubplot:xlabel='Place'>
К кластеру 2 типа относятся города из разных штатов в равной степени, за исключением двух: Флорида и Техас. Штаты, города которых входят в 0 кластер, также очень разнообразны, нельзя выделить какой-нибудь доминирующий, однако немного выделяется Техас. В отличие от двух предыдущих кластеров, во 1 кластере очевидно доминирующая доля Калифорнии, города которой занимают примерно 25% от общего количества. За ней следуют Коннектикут и Нью-Джерси, количество городов которых совокупно также составляет около 25%.
В целом, в 0 кластер входит большое количество штатов, что логично, ведь 0 кластер отражает "средние" города. Наоборот, в кластер 1 входит меньшее количество штатов.
m = folium.Map(zoom_start=50)
df['Long'] = data['Long']
df['Lat'] = data['Lat']
for i in data.index:
if df.loc[i, 'class'] == 0:
folium.Marker(
location=[df.loc[i, 'Lat'], df.loc[i, 'Long']],
icon=folium.Icon(color="blue"),
).add_to(m)
elif df.loc[i, 'class'] == 1:
folium.Marker(
location=[df.loc[i, 'Lat'], df.loc[i, 'Long']],
icon=folium.Icon(color="green"),
).add_to(m)
else:
folium.Marker(
location=[df.loc[i, 'Lat'], df.loc[i, 'Long']],
icon=folium.Icon(color="red"),
).add_to(m)
m
Из карты видно, что крупные города сосредоточены на западном побережье в Калифорнии (где находится большее количество айти-предприятий и развиты промышленность и киноиндустрия) и на восточном побережье вокруг Нью-Йорка -- финансовой столицы США.